

*This file conducts the estimation for the rho sensitivity (Table B.5)

*Bring in regression data
use regression_data_sensitivity_rho.dta, clear

*drop ROW as an originId
drop if originId==62

egen clusterId=group(originId)

************************************************************************************

*Collapsed version: average across importers
collapse logX_theta logL* instrument_* clusterId theta se_gyy, by(year sectorId originId) 
mkmat se_gyy if originId==1 & year==2010 & sectorId>2, matrix(S)
mkmat theta if originId==1 & year==2010 & sectorId>2, matrix(T)
**************************************************************************************
*Multiple Elasticities: Heterogeneous TE
**************************************************************************************

*FE
local dvars logX_theta
local indvars c.logL#sectorId

foreach x in low high cd	{

	*CES Instrument
	local instruments c.instrument_`x'#sectorId
	display "`x'"
	ivreghdfe `dvars'  (`indvars'=`instruments') if sectorId>2 & year==2010, cluster(clusterId) absorb(sectorId#year originId#year)
	*outreg2  using ${tab_dir}/donut.xml, append excel nor2 dec(2)  addstat(F-stat, e(widstat)) noaster nopa
	*Save statistics for table
	matrix A`x'=e(V)
	matrix b`x'=e(b)
	}

	
*Need to compute standard errors
keep if year==2010 & sectorId>2
keep originId sectorId logX_theta logL instrument_low instrument_high theta instrument_cd
sort sectorId

*First step: make matrix of instruments and independent variables, plus the matrix of the derivative of the dependent
*variable with respect to the trade elasticity vector

forvalues x=3/17	{
	gen logX_theta_`x'=-1/theta*logX_theta if sectorId==`x'
	replace logX_theta_`x'=0 if logX_theta_`x'==.
	gen logL_`x'=logL if sectorId==`x'
	replace logL_`x'=0 if logL_`x'==.
	gen instrument_low_`x'=instrument_low if sectorId==`x'
	replace instrument_low_`x'=0 if instrument_low_`x'==.	
		gen instrument_high_`x'=instrument_high if sectorId==`x'
	replace instrument_high_`x'=0 if instrument_high_`x'==.	
		gen instrument_cd_`x'=instrument_cd if sectorId==`x'
	replace instrument_cd_`x'=0 if instrument_cd_`x'==.	
	}

*Second step: residualize the logX_theta`x' and logL_`x' with respect to the fixed effects

forvalues x=3/17	{	
	reghdfe logX_theta_`x' if sectorId>2,  ///
	absorb(sectorId originId) res(logX_thetahat_`x') 
	reghdfe logL_`x' if sectorId>2,  ///
	absorb(sectorId originId) res(logLhat_`x') 
	
}

reghdfe logX_theta if sectorId>2,  ///
absorb(sectorId originId) res(y) 


*Third step: construct matrices of instruments and residualized independent and dependent variables
mkmat instrument_low_*, matrix(ZL)
mkmat instrument_high_*, matrix(ZH)
mkmat instrument_cd_*, matrix(ZCD)
mkmat logX_thetahat_*, matrix (Y)
mkmat logLhat_*, matrix(L)
mkmat y, matrix(y)

*Fourth Step: Compute the Jacobian matrices for OLS and IV and the reduced form
matrix JL = -inv(ZL'*L)*ZL'*Y
matrix JH = -inv(ZH'*L)*ZH'*Y
matrix JCD = -inv(ZCD'*L)*ZCD'*Y

*Compute the variance-covariance matrix and the standard errors

matrix VDL = vecdiag(Alow+JL*diag(S)*diag(S)*JL')
mata : st_matrix("SEL",  st_matrix("VDL"):^(.5))

matrix VDH = vecdiag(Ahigh+JH*diag(S)*diag(S)*JH')
mata : st_matrix("SEH",  st_matrix("VDH"):^(.5))

matrix VDCD = vecdiag(Acd+JCD*diag(S)*diag(S)*JCD')
mata : st_matrix("SECD",  st_matrix("VDCD"):^(.5))

*Preparing for export

matrix M = blow\SEL
mata: st_matrix("m1", colshape(st_matrix("M")', 30))
matrix M = bhigh\SEH
mata: st_matrix("m2", colshape(st_matrix("M")', 30))
matrix M = bcd\SECD
mata: st_matrix("m3", colshape(st_matrix("M")', 30))


matrix M=(m1\m2\m3)'
matrix list M

putexcel set table_rhonut, sheet(SE Estimates) replace
putexcel A1=matrix(M)
